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1*3  abstract  ] 

We  consider  the  usual  univariate  linear  model  E(y)  =  X7  ,  V(y)  =  o^I  . 

In  Part  One  of  tnis  paper  X  has  full  column  rank.  Numerically  stable  and 
efficient  computational  procedures  are  developed  for  the  least  squares  estima¬ 
tion  of  7  and  the  error  sum  of  squares.  We  employ  an  orthogonal  triangular 
decomposition  of  X  using  'Householder  transformations.  A  lower  bound  for  the 
condition  number  of  X  is  immediately  obtained  from  this  decomposition. 
Similar  computational  procedures  are  presented  for  the  usual  F-test  of  the 
general  linear  hypothesis  L'7  =  C  :  L'7  =  m  is  also  considered  for  m  ^  0  . 
Updating  techniques  are  given  for  adding  to  or  removing  from  (X,y)  a  row,  a 
set  of  rows  or  a  column. 

In  Part  Two,  X  has  less  than  full  rank.  Least  squares  estimates  are 
obtained  using  generalized  inverses.  The  function  L’7  is  estimable  when¬ 
ever  it  admits  an  unbiased  estimator  linear  in  y  .  We  show  how  to  computa¬ 
tionally  verify  estimability  of  L'7  and  the  equivalent  testability  of 
L'7  =  0  . 
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by 
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and 
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McGill  University 

Abstract 

2 

We  consider  the  usual  univariate  linear  model  E(y)  =  X7  ,  V(y)  =  a  I  . 
In  Part  One  of  this  paper  X  has  full  column  rank.  Numerically  stable 
and  efficient  computational  procedures  are  developed  for  the  least  squares 
estimation  of  7  and  the  error  sum  of  squares.  We  employ  an  orthogonal 
triangular  decomposition  of  X  using  Householder  transformations.  A  lower 
bound  for  the  condition  number  of  X  is  immediately  obtained  from  this 
decomposition.  Similar  computational  procedures  are  presented  for  the 
usual  F-test  of  the  general  linear  hypothesis  Lr7  =  0  ;  L'7  =  m  is 

also  considered  for  m  /  0  .  Updating  techniques  are  given  for  adding  to 
or  removing  from  (X,y)  a  row,  a  set  of  rows  or  a  column. 

In  Part  Two,  X  has  less  than  full  rank.  Least  squares  estimates  are 
obtained  using  generalized  inverses.  The  function  L'7  is  estimable 
whenever  it  admits  an  unbiased  estimator  linear  in  y  .  We  show  how  to 
computational] y  verify  est inability  of  L*7  and  the  equivalent  testability 
of  L*7  =  0  . 
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PART  ONE:  UNIVARIATE  LINEAR  MODEL  WITH  RJLL  RANK 

1.  Least  squares  estimr  'on  and  error  sum  of  squares 
We  consider  the  univariate  general  linear  model 

(1.1)  E(y)  =  X  7  5  V(y)  =  akI  , 

where  EC*)  denotes  mathematical  expectation  and  V(  •)  the  variance- 
covariar.ce  matrix.  We  take  the  design  matrix  X  to  be  nxq  of  rank 
q  <  n  and  known;  in  part  two  we  relax  this  assumption  of  full  column 
rank.  The  unknown  vector  y  of  q  regression  coefficients  is  estimated 
by  least  squares  from  an  observation  y  by  minimizing  the  sum  of  squares 

(1.2)  (y  -  X/)*(y  -  Xy)  . 

Prime  denotes  transposition;  bold-face  capital  letters  denote  matrices 
and  bold  lower-case  letters  vectors,  with  rows  always  appearing  primed. 

In  the  case  where  V(y)  =  a  A  in  (1.1),  with  A  known  and  positive 
definite,  we  may  replace  y  by  Fy  and  X  by  EX  where  F  satisfies 
FAF*  =  I  .  The  matrix  F  is  not  unique  but  it  is  possible  to  find  an  F 
which  is  lower  triangular  from  the  Cholesky  decomposition  of  A  (cf.  e.g., 
Healy,  1968) . 

It  is  well  known  that  the  least  squares  estimate  7  satisfies  the 
normal  equations 

(1.5)  X*X7  =  X’y 

and  is  unique  when  X  has  full  rank.  The  matrix  X*X  is  greatly 
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influenced  by  roundoff  errors  and  is  often  ill-conditioned :  by  this 

we  mean  that  a  relatively  "small"  change  in  X  will  induce  a  correspondingly 

"large"  change  in  (X^)"1  and  in  the  solution  y  =  (X’X)-1X*y  to  (1.3). 

For  these  reasons  we  prefer  to  work  with  X  directly  rather  than  X'X 
[cf.  e.g.,  Longley  (1967),  Wampler  (1969,  1970)]. 

It  is  possible  to  find  an  n  x  n  orthogonal  matrix  P  such  that 


(l.h) 


X  = 


P'X  = 


where  R  is  upper  triangular  of  order  q  x  q  .  This  ortho  ^nal  triangular 
decompos it ion  (CTD)  may  be  made  in  various  ways;  a  very  stable  numerical 
procedure  (Golub,  1965)  is  to  obtain  P  as  the  product  of  q  Householder 
transformations . 

A  square  matrix  of  the  form  H  =  I  -  2uu*  ,  where  u*u  =  1  ,  is  defined 

to  be  a  Householder  transformation.  Clearly  H  =  H*  and 
2 

HH’  =  H*H  =  H  =  I  ,  so  that  H  is  a  symmetric  and  orthogonal  matrix. 

All  but  one  of  the  characteristic  roots  of  H  are  unity,  the  simple 
root  being  -1  . 

A  vector  x  may  be  transformed  by  a  Householder  transformation  to 
a  vector  with  each  element  zero  except  for  the  first,  i.e.. 


(3-5)  Hjk  =  rei  ;  r  ^  0  , 


say,  where  e.  is  an  1* 1  vector  with  each  component  0  except  for 

~  <J 

the  j-th  which  is  1  f.i  =  l,2,...,n)  .  Premultiplying  (1.5)  by  its 
transpose  yields 

(1.6)  x’x  =  x'H'Hx  =  r2eje1  =  r2  . 
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Substituting  H  =  I  -  2uu*  in  (1.5)  gives 

(1.7)  x-2(u'x)u  =  re1  ; 

premultiplication  by  u*  yields  -ufx  =  ru..  ,  where  u  =  e’u  ,  the 

MW  J.  X  wlw 

first  element  in  u  .  Substitution  in  (1.7)  gives  x+2ru^u  =  re^  , 

so  that  with  x  =  {x^}  , 

(1.8)  2u2  =  l-(x1/r)  ;  2u±  =  -x^ru^)  ,  i  =  2,...,n  . 

The  first  expression  will  always  be  computed  positive  if  the  square  root 
of  (1.6)  is  taken  as 

(1.9)  r  =  -sgn(x1) -(x'x)1/1-  , 

where  sgn(x1)  =  +1  if  x.^  >  0  and  -1  otherwise.  Then 

(1.10)  2ux  =  1+  ( Ix-J/s)  ;  2ui  =  sgn(x1)*xi/(su1)  ,  i  =  2,  ...,n, 

where 

(1.11)  s  =  +(x*x)1//2  . 

This  gives  u»u  =  1  ,  f or  2u2  =  x2/ (2s2u2)  =  x2/(s2+  s  |x1|)  ,  i  =  2,  ...,n  . 

Hence  2£u2  =  (s2"xi)/(s2+  s  lxql)  =  1"(|X1I/B)  =  2(l-u2)  .  We  note 
i=2 

that  H  need  not  be  confuted  explicitly  as  Hx  =  x  -  2(u’x)u  ,  for  which 

we  need  only  u  and  v'  -  .  In  the  above  form,  it  is  necessary  to  compute 

two  square  roots  per  Householder  transformation;  if,  however,  we  write 
H  =  I  -u(u,u)’1u'  then  only  one  square  root  need  be  calculated  (Businger 
and  Golub,  1965) . 
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we  obtain 


Applying  this  procedure  with  x  replaced  by  Xen  , 

(1.12)  HX  =  (r^y  , 

where  r^  replaces  r  ,  and  is  nx  (q-1)  such  that 

X,e.  =  x.j,  -  2(u'x.^.)u  ,  j  =  l,...,q-l  and  x.±1  =  Xe.J_  .  This 
~l~j  -j+1  x _ j+1  „  ~j+l  ~~j+l 

procedure  is  now  repeated  with  as  x  and  a  Householder  transformation 

H1  =  I  -  2u^u^  ,  say,  with  =  0  .  The  last  n-2  elements  of  X^e^ 

are  now  annihilated.  So  H-jHXe.^  =  r11H.|e.|  =  rnei  3  while  H^HXeg  =  Hl^ifl 
has  its  last  n-2  components  zero.  The  product  H^H  is  orthogonal. 

Further  repetitions,  annihilating  at  the  j-th  stage  the  last  n-j  elements 
in  the  j-th  column  of  the  matrix  X  transformed  previously  by  j-1 
Householder  transformations  ( j  =  1,  ...,q)  ,  realizes  P  as  the  product 
of  q  Householder  transfoimations.  The  matrix  P  is  not  computed 
explicitly.  Details  of  this  algorithm  are  given  by  Golub  (1965),  and 
Businger  and  Golub  (1965)  who  also  give  a  program  in  Algol  60. 

Partitioning  P  =  (P^Pg)  ,  with  P^  nxq  and  Pg  nx(n-q)  gives 
from  (1.4) 

(1.13)  PjX  -  R  ,  F*X  =  0  , 


with  P»Pn  =  I  ,  P*P0  =  0  and  P»P  =  I  ,  since  P'P  =  I  .  If,  in 
<v Ik  1  «w(^  mXm ~n-q  «v  ^  <vH 

the  above  algorithm,  we  simultaneously  apply  the  q  Householder  transfor¬ 
mations  to  the  observation  vector  y  ,  then  we  have 


(l.l4)  P»y  = 


f  b 

1st  J  h 


say.  Thus  Zg  =  P£y  has  expectation  E(P^y)  =  P^X>  =  0  and  covariance 


4 


2  2 

matrix  V(P^y)  =  o  P^P2  =  a  ln_  .  Hence  z 2  is  an  easily  computed 
vector  of  uncorrelated  regression  residuals  and  may  be  used  to  test  for 
serial  correlation  (cf.  e.g.,  Grossman  and  Styan,  1970).  It  follows  that 

(1.15)  +  x(x‘x)_1x«  =  in  , 


as  each  term  on  the  left-hand  side  is  idempotent  and  their  cross-product 
is  0  ;  their  sum  is  idempotent  with  rank  the  sum  of  the  ranks  n-q  and  q 
So  PgP*  =  I  -X^y)"1^  and  2«zp  =  y*P?P«y  =  (y  -  Xy)  » (y  -  Xy)  is 
the  error  sum  of  squares  3g  ,  say  —  the  minimum  of  (1.2)  .  It  is  simply 
computed  here  as  the  sum  of  squares  of  the  n-q  elements  in  z2  =  P^y  . 

A 

The  vector  of  (correlated)  residuals  r  =  y-X7  is  often  essential 
for  analysis  of  the  linear  model  (cf.  e.g.,  Draper  and  Smith,  1966) . 

Though  the  matrix  P  may  not  be  computed  explicitly  it-  can  be  retrieved 
as  the  product  of  q  Householder  transformations  when  the  corresponding 
q  u  vectors  have  been  stored  (which  we  recommend) .  Hence  we  compute 
r  =  P2z2  ,  since  ?2z2  =  PgP^y  =  [I  - X(X*X) ""Sc*  ]y  =  y-X7  =  r  .  However, 
it  has  been  observed  by  Gentleman  (1970)  that  computing  r  in  this  fashion 
may  be  numerically  unstable. 

We  also  find  from  (1.4)  that 


(1.16)  X»X  =  (R 


=  R'R  . 


Substitution  in  (1.5)  yielus  R'R7  =  (R,,0)P,y  =  R’z^  >  so  ^4at  solving 


(1.17) 


R7  =  z. 


gives  7  .  This  is  expedited  by  R  being  upper  triangular. 


We  note  that  R*R  is  a  Cholesky  factorization  of  X*X  ,  for  which 

Healy  (1968)  has  given  a  Fortran  program. 

*  -2-1 
The  estimator  7  has  covariance  matrix  V(7)  =  a  (X*X)  ;  an 

unbiased  estimate  is  Sg(X,X) "^/(n-q)  which  is  easily  computed  using 

(1.16)  as  (z^z2)r“1(R-1) '/(n-q)  .  The  generalized  variance  (cf.  e.g., 

Anderson,  1958)  is  |v(7)  j  =  a^q/|x’X|  ,  where  {•  J  denotes  determinant. 

In  optimal  design  theory  a  problem  is  to  choose  X  so  that  |X’X  |  is 

maximized  thus  reducing  |v(7)  |  as  much  as  possible.  Again  using  (l.l6) 

q  2 

we  see  that  |X’X|  =  |r*R|  =  |  |  rv.  ,  as  R  is  upper  triangular.  Hence 

i=l 

|V(7)  |  is  estimated  by  [ z£z0  /  (n-q)  ]q  /  "TV  a  • 

i=l 

A  measure  of  the  ill-conditioning  of  a  matrix  is  its  condition  number 
which  we  define  as  the  ratio  of  the  largest  and  smallest  nonzero  singular 
values  of  the  matrix.  The  singular  values  of  a  (possibly  rectangular) 
matrix  A  are  the  positive  square  roots  of  the  characteristic  roots  of 
A*A  or  AA*  .  When  the  condition  number  far  exceeds  the  rank  we  find 
(cf.  Wilkinson,  1967)  that  the  matrix  is  extremely  ill-conditioned. 

A  lower  bound  for  the  condition  number  x(X)  of  the  design  matrix  X 
is  the  ratio  of  the  largest  and  smallest  (in  absolute  value)  diagonal 
elements  of  R  .  To  see  this  we  note  first  that  X  and  P'X  have  the 
same  singular  values,  due  to  the  orthogonality  of  P  .  As  P'X  is  merely 
R  bordered  by  zeroes,  sg(X)  =  sg(R)  ,  where  sg(*)  denotes  singular 
value.  For  any  square  matrix  A  of  order  n  x  n  , 

(1.18)  sgn(A)  <  |chj(A)|  <  sg^A)  ;  i  =  1,  ...,n 

with  ch(*)  denoting  characteristic  root.  The  subscript  j  indicates 
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O-th  largest.  To  prove  (1.18)  when  A  has  real  roots,  let  \  =  ch.(A) 

<3  ~ 

with  Av  =  \v  .  Then 
2 

(1.19)  sg1(A)  =  ch-^A’A)  =  max[x'A,Ax/x,x]  >  v’A’Av/v’v  = 

x  ~  ~  ~  ~  ~~  ~  ~ 

=  \v’A,v/v,v  =  X2 
Similarly  sg2(A)  <  \2  .  Thus 


sg-.  (X)  sg  (R)  maxlch(R)  I  max|r..| 

(1.20)  k(X)  =  >  min  |ch(R)|  =  SinjTTTJ  * 

Other  properties  of  k(A)  are  given  by  Wilkinson  (1967). 

Why  is  the  condition  number  important  and  how  can  we  use  the 
relationship  (1.20)?  Let  7  be  the  computed  approximation  to  7  which 
satisfies  (1.3) •  Suppose  that  we  wish  to  determine  an  upper  bound  for 
the  norm  of  the  relative  error  of  7  : 

(1.21)  |!  7  -  7  1| /||  7  ||  , 

where  ||  a  ||  indicates  the  Euclidean  norm  (a’a)1^  .  Define 

(1.22)  r  =  y  -  Xy  , 


which  we  can  compute  quite  accurately.  Then 


(1.23)  r  -  r 


y 


and  hence 

(1.24)  -X’X(  7-7)  =  -X»?  , 


7 


since  X*r  =  0  .  Thus 

(1.25)  ||7 -7  ||  =  ||(X-X)"1X'r||  <  ch1[(X*X)‘1]||X'r||  =  ||X'r||/sg^(X)  . 

From  (1.3),  1|X'X7||  =  ||X'yi|  ,  so  that 

(1.26)  |jX'yl|  <  ||  7  ||sg^(X)  . 

Combining  (1.25)  and.  (1*26),  we  have 

(1.27)  ||  7 -7  ||  /  ||  7  ||  <  f sg!^)  /  sSq[  (X)  ]2  ||  X»r  ||  /  ||  X»y  || 

=  K2Wl|X,r||/||X,y||  . 

Thus  we  see  that  the  condition  number  may  be  used  for  determining  an 

upper  bound  for  the  relative  error  of  ||  y  ||  .  This  upper  bound  is  the 

2 

product  of  two  factors;  the  first  of  which,  k  (X)  ,  is  independent  of  y 
However,  the  lower  bound  provided  by  (1.20)  would  in  some  circumstances 
give  insight  into  the  relative  error.  Hence,  if 

(1.28)  [max  |rii  |  /  min|rijL  |  f  |{  X'r  ||  /  ||  X’y  || 

is  large,  then  it  is  likely  that  the  relative  error  in  ||7  ||  is  large. 

The  numerical  efficiency  of  the  above  orthogonal  triangular 
decomposition  is  enhanced  (cf .  Golub,  1965)  if  the  column  selected  for 
each  of  the  q  Householder  transformations  maximizes  the  corresponding 
sum  of  squares.  That  is,  at  the  j-th  stage  (j  =  1, ...,q)  we  transform 
that  column  of  the  q-j+1  possibilities  which  maximizes  the  sum  of 
squares  of  its  last  n-j+1  components.  The  interchanges  may  be 
summarized  in  a  permutation  matrix  Tf  postmult iplying  X  .  Thus  (1.4) 


becomes 


The  vector  z  does  not  change  and  hence  neither  does  Sg  .  The  solution 
(1.17)  changes  however;  substituting  (I.29)  into  (1.3)  now  gives 
TTR*RTr'7  =  TTR'z1  ,  so  that 


(1.30)  R(TT’7)  =  z±  =  R9  , 


is  solved  for  9  ,  and  7  =  TT9  .  As  these  interchanges  only  rearrange 

q  2 

the  r..  we  still  find  Ix'xl  =  TTr..  .  The  lower  bound  for  the  condition 

xi  ' _ 1  !  ’  11 

i=l 

number  simplifies,  however,  as  with  these  interchanges  max  |r^  |  =  |r^  |  , 
and  min|ri;.|  =  |rqq|  so  that  k(X)  >  |rn/rq<ll  • 

Given  the  n  x  n  matrix 

1  ,  "1  ,  “1  ,  ...  ,  — 1 
0,  1,-1,. ..,-1 
(1.31)  a  =  :  :  :  : 

~  f  •  •  • 

0  y  0  j  0  j  j  1 


we  see  that  max|r.„  j  =  min|r._  |  =  1  ,  and  so  k(A)  >  1  ,  since  A  =  R 
when  no  column  interchanges  are  made.  However,  if  column  interchanges 
are  performed  then  for  u  -  10  say,  |r^j  =  3*162,  jr^j  A  .003383 
and  k(A)  >  93^.8.  The  actual  value  of  k(A)  =  1918.5  • 

The  Fortran  IV  programs  LLSQ  and  DLLSQ,  (double-precision)  in  the 
Scientific  Subroutine  Package  (SSP)  of  IBM  (1968)  solve  the  least  squares 
problem  as  described  above.  The  SSP  library  is  available  at  many  IBM  360 
computing  centers.  The  SSP  manual  gives  a  write-up  of  the  procedure  and 
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indicates  how  y  and  Sg  are  output.  In  addition  we  note  that  the  q 
diagonal  elements  of  R  are  output  as  '  AUX(q+l, . .  .,2q)  *  ,  with 
raaxlr.^1  =  AUX(q+l)  and  min|r.„J  =  AUX(2q)  in  absolute  value.  The 
remaining  nonzero  elements  of  R  are  overwritten  in  corresponding 
positions  of  X  (input  as  *  A  ’)•  The  vector  z  is  overwritten  on  y 
(input  as  *  B  ')  and  Sg  appears  in  *  AUX(l)  ’ .  The  solution  y  is 
output  as  f  X  ' . 

2  3 

The  number  of  multiplications  to  obtain  R  is  about  nq  -  q"/3  , 

2 

whereas  approximately  nq  /2  multiplications  are  required  to  form  the 
normal  equations  (1.3)  with  about  q^/6  multiplications  needed  to  solve 
them.  Thus  when  n-q  is  small,  the  number  of  operations  is  roughly  the 
same  for  both  algorithms,  but  when  n-q  is  large,  it  requires  about  twice 
as  many  operations  to  use  the  orthogonalization  procedure. 

The  orthogonal  triangular  decomposition  (1.4)  or  (I.29)  is  very 
similar  to  the  Gram-Schmidt  decomposition.  Indeed  if  n  =  q  and  there 
is  no  roundoff  error  and  all  r^  are  taken  positive,  then  the  Householder 
and  Gram-Schmidt  algorithms  yield  precisely  the  same  transformation. 
Although  the  modified  Gram-Schmidt  process  (cf.  e.g.,  Golub,  1969)  may  be 
used  for  solving  linear  least  squares  problems,  the  computed  vectors  may 
not  be  truly  orthogonal!  The  Householder  transformations,  however,  yield 
vectors  which  are  more  nearly  orthogonal  (Wilkinson,  1965) .  Furthermore, 
not  only  do  the  first  q  columns  of  P  span  the  same  space  as  the 
columns  of  X  ,  but  the  last  n-q  columns  of  P  span  the  complement  of 
the  space  spanned  by  the  columns  of  X  .  As  we  have  seen  above,  this  is 
quite  useful. 
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Hypothesis  testing  and  estimation  under  constraints 
Let  us  consider  the  general  linear  hypothesis 

(2.1)  L'r  =  0 

for  the  linear  model  of  Section  1.  The  contrast  matrix  L*  is  taken  as 
s  x  q  of  full  row  rank  s  <  q  .  If  we  assume  that  y  is  normally- 
distributed  then  L’7  is  N(L’7>  c/Tj*  (X'X)’"^)  ,  with  y  =  (X,X)~'S{,y  . 
The  numerator  of  the  usual  F-test  for  (2.1)  is  then  well  known  to  be 

(2.2)  7 rL[Lr (X'X) _1  L  J”1!'  7  =  Sh  , 

say,  the  "hypothesis  sum  of  squares".  Substituting  (1.16)  and  (1.17) 
into  (2.2)  gives 

(2.3)  Sh  •  z^R'VLtL'R’VVL]'1!*  R-^  . 

We  compute  (R  ^) 'L  =  G  ,  say,  by  solving  R’G  =  L  ,  with  Rf  lower 
triangular.  We  then  obtain  an  orthogonal  triangular  decomposition  of  G  , 
qx  s  (q  >  s)  , 

(2.4)  G  =  (R"1)»L=Q 

say,  where  B  is  upper  triangular  s  x  s  and  the  orthogonal  matrix  Q 
is  the  product  of  s  Householder  transformations.  Then  G'G  =  B’B  ; 
partitioning  Q  =  (Q-^Qg)  ,  where  Q1  is  qxs  and  Qg  qx  (q-s)  gives 
G  =  Q^B  from  (2.4).  Substitution  in  (2.3)  yields 

(2.5)  Sh  = 
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which  we  compute  by  applying  the  s  Householder  transformations  of  (2.4) 
to  simultaneously  with  G  and  then  summing  the  squares  of  the 

first  s  components  of  the  transformed  . 

If  we  test  the  hypothesis 

(2.6)  L'y  =  m  , 

where  m  is  a  given  s  x  1  vector,  not  necessarily  0  ,  then  we  proceed 
by  computing  L xy  -  u  =  h  ,  say,  and  sum  the  squares  of  the  components  of 
(B  ^)’h  ;  we  find  the  latter  by  solving  L *y  -m  =  h  =  B*t  ,  say,  for  t  , 
with  B*  lower  triangular. 

The  described  procedure  can  be  improved  upon  when  s  >  q-s  .  We 
first  obtain  an  orthogonal  triangular  decomposition  of  L  , 

(2.7)  L  =  T 

say,  where  T  is  orthogonal  and  U  upper  triangular.  Partitioning 
T  =  (T^Tg)  ,  where  T  is  qxs  and  T2  is  qx  (q-s)  leads  to 

(2.8)  L,T1  =  U’  ;  L'Tg  =  0  . 

Thus  L*7  =  0  if  and  only  if  7  =  Tg©  for  seme  9  ,  now  unconstrained. 
Hence 

(2.9)  min  (y-X7)*(y-X7)  =  min(y -XTp9)*(y -XT„9)  . 

L»7=0  ~  ~~  "  ~  9  ~  ~~~  ~ 

Using  (1.4)  and  (1.14),  we  see  that  (2.9)  reduces  to 

(2. ID)  m^^-HT^^^-RTg©)  +  z^z2  , 

0 
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so  that  equals  the  first  term  in  (2.10)  which  is  easily  computed 

as  in  Section  1  with  replacing  y  and  RT2  replacing  X  .  Since 
(cf.  e.g.,  Good  (1965),  p.  89)> 


(2.11)  sg  (XT)  >  sg  (XT  )  , 

1  1  MMW 


(2.12)  sg^XT)  <  sgq.s+1(XT)  <  seq.s(®2)  ' 

we  have 

(2.13)  k(XT2)  <  h (XT)  =  h(X)  . 

Thus,  by  eliminating  the  constraints,  the  linear  least  squares  problem 
may  become  better  conditioned. 

* 

The  least  squares  estimate  7  ,  say,  of  7  subject  to  L*7  =0  is 

obtained  from  the  solution  0  to  (2.10)  by 

(2.14)  7*  =  T29  . 


If  the  constraints  have  nonnull  righthand  side  m  as  in  (2.9)  then 
the  procedure  is  changed  as  follows.  Evidently  1/7  =  m  holds  if  and 
only  if  7  =  TgQ  +  T^U  'L)  *m  =  TgQ  +  T^  ,  say.  We  obtain  w  by  solving 
m  =  U’w  ,  with  U*  lower  triangular.  Thus  y  is  replaced  by  y  -  XT^w 
and  hence  z1  by  z^  -  RT^w  the  resulting  value  of  is  therefore 


(2.15)  min(z  -RT  w  -  ,,T  0)»(z  -RT  w-RT~©) 


which  we  compute  as  in  Section  1  with  z^  -  RT^w  replacing  y  and  RTg 


replacing  X 


The  relevant  F-test  for  the  hypotheses  (2.1)  or  (2.6)  is  then 
computed  as 

Vs 

(2'16)  F  -  WJWV  ■ 

with  the  critical,  region  formed  by  values  of  (2.l6)  exceeding  the  corres¬ 
ponding  tabulated  value  of  F  with  s  and  n-q  degrees  of  freedom. 

In  some  special,  though  common,  situations  the  above  computations 
simplify  considerably. 

If  we  test  a  single  contrast  in  7  equal  to  0  we  obtain  (2.1) 
with  s  =  1  .  Let  us  write  this  as 

(2.17)  1*7  =  0  . 

A  particular  case  might  be  testing  a  single  regression  coefficient  equal 
to  0  .  Then  (R_1) *L  =  K  becomes  (R-1)*!  =  k  ,  say,  found  by  solving 
l  =  R'k  as  before.  Then  (2.3)  becomes 

(2.18)  (l,7)2/k,k  =  sh  , 

and  we  compute  the  denominator  in  (2.l8)  by  summing  squares  of  components 
in  k  .  The  one-sided  t-test  for 

(2.19)  1*7  >  0 

has  critical  region  large  positive  values  of  I* 7/  [krk  S  /(n-q) l1/2  . 

Another  special  case  occurs  with  s  =  q-1  when  L'7  =0  if  and 
only  if 

(2.20)  7  =  9t  , 
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where  9  is  now  a  scalar.  The  vector  t  is  often  found  upon  inspection 
(without  transforming  L  ) .  For  example  in  testing  for  homogeneity  of 
coefficients  of  y  ,  we  have  t  =  e  ,  the  vector  with  each  component 
unity.  Substituting  t  for  in  (2.10)  yields 

(2.21)  9  =  z^Rt /t»R*Rt  , 
and 

(2.22)  Sh  =  "  (z-jRt)2  /  t'R'Rt  , 

with  the  denominators  computed  by  summing  squares  of  elements  of  Rt  . 
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After  a  particular  set  cf  data  has  been  analyzed  it  is  often 
pertinent  to  add  to  or  remove  from  X  and  y  a  row  (or  set  of  rows 
or  to  add  to  or  remove  from  X  a  column.  This  happens  when  new  informa¬ 
tion  becomes  available  or  when  existing  experimental  units  have  been 
classified  as  extreme,  or  independent  variables  insignificant. 

We  begin  by  considering  the  addition  of  data  from  m  ,  say,  further 
experimental  units.  Let  X  and  y  be  the  corresponding  data  of  order 

rJTl  rJU 

mxq  and  mxl  respectively.  Following  (1.4)  and  (1.14)  we  may  write 


(3.1) 


f  I 

0  ^ 

f  x 

y 

~m 

lm 

0 

p* 

X 

y  , 

-  J 

V  ~ 

~  J 

Applying  q  Householder  transformations  of  order  ntf-q  to  the  first  m+q 
rows  of  (3.1)  yields 


-ft 

say,  where  is  qxq  upper  triangular,  zQ  is  q x  1  and  z 1  is 
mxl.  Hence 


1 6 


(3.U) 


0 


jn 


I 

_n-q 


J 


is  an  orthogonal  matrix  formed  from  2q  Householder  transformations,  and 

has  order  ntf-n  .  The  new  residual  sum  of  squares  is  z*’z*  +  Z2Z2  * 

i.e.,  the  previous  sum  of  squares,  z^z2  ,  augmented  by  the  sum  of  squares 

* 

of  the  m  component s  of  z1  ;  these  components  themselves  give  m 
additional  uncorrelated  residuals. 

Next,  suppose  we  wish  to  add  a  (q+l)-tb  variable  whose  n  values 
constitute  a  vector  x  .  We  first  compute  P’x  by  applying  in  turn  the 
q  Householder  transformations  determined  by  the  stored  vectors  u 
(cf.  residual  calculations  in  Section  1).  We  need  then  only  one  further 
Householder  transformation,  H  ,  say,  of  order  n-q  to  annihilate  the 
last  n-q-1  elements  in  P’x  ,  i.e.. 


rh 

f  5  Hi  ) 

f  R  PjxA 

P’(X,x)  = 

= 

0  H  I 

1  0  HF'x  J 

l  0  he..  , 

~  J 

'v.  -  ~2"  J 

y.  ~  -1-/ 

where  P  =  (P^Pg)  ,  as  in  §1,  and  h  =  x’PgP^x  ,  the  sum  of  squares  of 
the  last  n-q  components  of  P’x  . 

The  procedure  for  removing  an  experimental  unit  is  more  complicated. 
The  method  given  previously  by  Golub  and  Saunders  (1970),  may  under 
certain  circumstances  prove  unstable.  We  now  give  a  new  method  which 
should  provide  a  more  accurate  solution. 
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Suppose  we  want  to  remove  x*  ,  the  i-th  row  of  X  .  We  seek  an 
upper  triangular  matrix  S  ,  say,  so  that 

(3.6)  x*x  "xi^i  =  R,R -x±x[  ■  s’s  =  RI(I -tt»)B  > 

say,  where  R*t  =  xi  ;  the  vector  t  is  easily  canputed  since  R*  is 
lower  triangular.  We  now  construct  an  orthogonal  matrix  Q  so  that 
Qt  =  ce^  ;  thus  c2  =  t*t  =  x! (R*R)~^x^  =  e’X(X,X)~1X,ei  <  1  .  We  define 
the  quasi-diagonal  matrices  of  order  q  x  q  : 

k  =  1, . .  »,q-l  , 


k  —  •  *  .,q-l  • 

and  9^  are  orthogonal.  Let 

(3*9)  =  Zc  -t  ^ -1,> -1^  *  *  =  t 

with  t_  =  t  and  Rrt  =  _  .  We  choose  9,  so  that  Z  .  annihilates 

~0  ~0  k  ~q -l 

e1  ,^-,t,  ..  and  hence  e*  ,,_t.  =  0  ;  l  =  1, ...,q-l  .  Then  the  matrix 
„q-!Ti_l-l  „q-z+l_z 

(5.1C)  S-5&—5,.! 

satisfies  Qt  =  ce^  and  is  orthogonal.  From  (3.6)  we  may  write 


Clearly 


(3.7)  Z*  « 


ik-1 


!k 


Jq-k-1 


where 


r 


(3.8)  ek  = 


cos  9,  ,  sin  9, 
k  7  k 


^  -sin  9k ,  cos  9k 
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(3*11)  S*S  =  R»Q‘(I -c2eie»)  QR  , 

p 

which  is  positive  definite  if  and  only  if  c  >  1  .  It  follows  that 


(3.12)  QR  =  W  = 


W11  '  W12  '  ***  Wl,q-1  '  Wlq 


W21  *  W22  *  *•*  W2,q-1  '  W2q 


0  , 


0  , 


W32  '  W3,q-1  '  W3- 


, . .  w  ,  ,  w 
q,q-l  qq 


is  an  upper  Hesseriberg  matrix.  Thus  (3.11)  becomes  S*S  =  W’D^  ,  with 
f (l-c2)1/2  0* 


(3.13)  D  • 


Jq-1 


which  is  real  when  c  <  1  .  We  compute  S  by  applying  orthogonal 
transformations  to  the  upper  Hessenberg  matrix  DW  .  Let 


(3.14) 


5k  = 


5k  5k-: 


k  —  1^  •  •  »f  q-1 


with  SQ  =  EW  and  formed  as  in  (3.7)  but  with  9^  replacing 

©k  and  so  chosen  that  z£  annihilates  £k  ek  =  ektlDWek  and  thus 

Jfcfl Ski  =  0  '  Then 

.  *  *  *  * 

(3.15)  s  =  S  x  =  Z  ±  Z^_2  ...  Z2  ZXDW 

p 

This  procedure  requires  about  9q  /2  multiplications  and  2q-l  square 
roots. 
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The  above  algorithm  can  also  be  used  for  adding  an  observation  but 
about  twice  as  many  numerical  operations  are  required  as  in  the  procedure 
given  by  (3.3)  and  (3  A)  •  We  also  note  that  the  problem  of  deleting  an 
observation  is  numerically  delicate.  Since 

(3.16)  S'S  =  R*(I  -tt»)R  , 
it  follows  that 

(3.17)  k(S)  <k(R)  /  (l-t*t)l/2  . 

Thus  if  t’t  is  close  to  1  ,  then  k(S)  could  be  quite  large  as  the 
right-hand  side  of  (3.17)  is  attainable. 

Finally  suppose  we  wish  to  remove  an  independent  variable  or 
column  of  X  .  If  it  is  the  last  then  no  further  calculations  are 
required;  but  suppose  it  is  the  first.  Let 


where  R  is  q  x (q-l)  and  has  one  more  row  than  an  upper  Hessenberg 
matrix.  We  annihilate  the  elements  just  below  the  main  diagonal  of  R  , 
i.e.,  t22}  rqq  '  by  pbplying  orthogonal  transformations  of  the  type 

(3.7)  with 


(3.19)  R*  =  ^  R^  5  k  =  l,...,q-l  , 
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and  R_  =  R  ;  we  choose  ©,  in 
_0  „  k 

is  annihilated;  thus  e*  n  R,  e, 

_k+l~k  „k 

matrix  sought. 


5k 


=  0 


so  that  e 
and  R 

~q-i 


k+1  ?k-l  fk  ~  rk+l,  k+1 
is  the  new  triangular 
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PART  TWO:  UNIVARIATE  LINEAR  MODEL  WITH  LESS 


THAN  FULL  RANK 

4.  Least  squares  estimation  and  error  sum  of  squares 

We  consider  now  the  univariate  general  linear  model  (1.1), 

(4.1)  E(y)  =  X7  ,  V(y)  =  a2 1  , 

with  the  design  matrix  X  of  rank  r  <  q  <  n  .  We  obtain  the  same  normal 
equations  as  (lo), 

(4.2)  X*X>  =  X’/ 

which  are  consistent;  their  solution,  however,  may  not  be  unique.  Consider 
a  solution  to  (4.2)  which  we  may  write 

(4.3)  7  =  (X’X)'X’y  , 

where  (•)  denotes  generalized  inverse.  We  .Follow  Pringle  and  Rayner  (1971) 

and  define  a  generalised  inverse  of  a  matrix  A  ,  m  x  n  ,  as  any  matrix  A 
satisfying 

(4.4)  AA~A  =  A 

Evidently  A  has  order  n  x  m  .  Such  a  generalized  inverse  exists  but  is 
not  unique  in  general;  i^;  however,  k  satisfies  (4.4)  and 

(4.5)  A" A  A"  =  A"  . 

(4.6)  ( AA  ) '  -  AA  } 

(4.7)  (A-AJ'.A-A  , 


y 


then  we  write  A  =  A  ,  the  pseudo -inverse  of  A  .  When  we  only  require 
that  (4. 4)  is  satisfied  we  will  write  A~  =  g^(A)  —  a  g^- inverse  of  A  . 

Similarly  when  (4.4)  and  (4.5)  are  satisfied,  A~  -  S-j^CA)  >  (4.4),  (4.5), 
and  (4.6):  A~  =  g^^{A)  .  The  pseudo-inverse  A+  =  g^^(A)  .  The 

A  A  A 

solution  ,  say,  to  (,L.2)  which  minimizes  7*7  equals  X  y  as  is 
shown,  for  example,  by  Peters  and  Wilkinson  (1970) •  Our  concern,  however, 
focuses  :.iore  on  estimable  functions  of  7  ,  rather  than  7  per  se  so  we 
will  not  discuss  here  computation  of  7Q  .  We  define  an  estimable  function 
of  7  as  a  vector  L'7  which  admits  an  unbiased  estimator  of  the  form 
K'y  ,  where  L‘  is  sxq  ,  say,  and  K*  ,  sxn  .  The  least  squares 
estimate  is  then  L*7  =  L’(X’X)*X,y  so  that  K*  =  L,(X*X)"Xr  .  We  shall 
see  (Section  5)  that  when  L'7  is  estimable,  L,(X*X)_X'  is  unique  for 
all  (X’X)  =  g^(X’X)  .  Rather  than  form  X'X  ,  find  a  g^X'X)  and  then 

postmult iply  it  by  X'  ,  we  compute  a  g12^(X)  directly,  noting  that  G 
is  a  g1(A)  if  and  only  if  it  can  be  written  as  (A'A)  A’  for  some 
g1(A*A)  =  (A’A)  (Pringle  and  Rayner  (1971) ,  P-  26], 

We  proceed  as  in  Section  1  to  orthogonally  transform  X  by  Householder 
transformations  with  column  interchanges.  If  X  has  rank  r  then  after  r 
Householder  transformations  we  obtain,  cf.  (I.29), 


R  S  \  /  R  S  \ 

(4.8)  X  =  p  ~  ~  I tT:  ;  P’  XTT  =  ~  ~  , 

~  ~  0  0  ~  ~  ~~  loo/ 


where  R  is  upper  triangular,  rxr  ,  S  is  rx  (q-r)  ,  and  tT  is  a 
permutation  matrix  of  order  qxq  .  We  now  claim  that 


(M) 


* 

X  =  Tr 


r"1  o\ 

~  P'  =  g,p,(x)  . 

0  0  ~  ^  ~ 
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7  =  X  y  to  (4.2)  afforded  by  (4.9)  is  often  called  a  basic  solution  as  it 
contains  at  most  q-r  nonzero  elements. 

Thus  (4.9)  accomodates  our  purposes;  moreover  we  do  not  have  a 
stronger  g-inverse  than  is  needed.  As  in  Section  1  we  partition 
P  =  (P1,P2)  ,  but  now  let  P1  be  nxr  and  Pg  nx(n-r)  .  From  (4.8) 
it  follows,  cf.  (1.13),  that 

(4.10)  PJXTT  =  (R,S) 


(4.11)  P»X  =  0 
Following  (1.14)  we  now  write 


where  z^  is  now  r  x  1  and  Zg  (n-r)  x  1  •  Thus  Zg  is  again  a  vector 
of  uncorrelated  residuals;  moreover 

(4.-13)  PgP£  +  £(£*£)  =  In  > 

as  in  (1.15),  with  Pg  P£  idempotent  rank  n-r  and  X^'X^X*  symmetric 
idempotent  rank  r  .  By  (4.11)  their  cross-product  is  0  and  so  their 
sum  is  idempotent  rank  (n-r)+r  =  n  and  hence  I  as  claimed.  Thus 
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(4.14) 


z*z2  =  y*(l-X(X*X)"X*)y 

is  the  residual  sum  of  squares,  computed  as  the  sum  of  squares  of  ti. 
n-r  components  in  • 

The  vector  of  (correlated)  residuals  r  =  y-X?  =  (I  -  X(X'X)  ~X')y  =  P?P£y 
as  in  Section  1,  and  using  (4.13)  it  follows  that  (4.l4)  equals  r'r  . 
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5.  Estimating  estimable  functions  and  testing  testable  hypotheses 
As  mentioned  in  Section  U  we  are  not  directly  concerned  with  the 
estimation  per  se  of  7  .  We  define  L'y  to  be  an  estimable  function 
of  7  whenever  it  admits  an  unbiased  estimator  which  is  linear  in  y  , 
K*y  ,  say.  Thus 


(5.1)  L'7  =  E(K’y)  =  K' X  7 


holds  for  all  y  .  Hence 

(5.2)  L*  =  K’X  . 


As  in  Section  3  we  take  L*  to  be  s  x  q  >  but  now  relax  the  assumption  of 
full  row  rank  taking  r(L)  =  t  <  r  .  We  obtain 


(5.3) 


9 


directly  from  (5*2).  Substituting  (4.8)  into  (5-3)  gives 


(5.10 


> 


where  we  partition 


(5-5)  L‘TT  =  (L*,L*)  , 

with  sxr  ,  and  s  x  (q-r)  •  The  matrix  L’TT  is  the  contrast 

matrix  L’  with  its  columns  permuted  according  to  the  interchanges  which 
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rearrange  the  columns  of  X  to  make  the  first  r  columns  linearly 
independent.  Then  L£  are  the  corresponding  r  columns  of  L'  or  L*TT  • 
We  apply  v  >  r  Householder  transformations  of  order  s+r  ,  whose 
product  is  V*  ,  say,  so  that 


where  TT^  is  a  permutation  matrix,  and  T  is  upper  triangular  vxv  . 

If  (5*6)  is  achieved  at  the  r-th  stage,  i.e.,  v  =  r  ,  then  L’y  is 
estimable.  If  not,  then  L*7  is  not  estimable. 

An  alternative  procedure  which  is  often  easy  to  verify  theoretically 
follows  and  is  included  for  completeness. 

THEOREM  $.1.  The  function  L'7  is  estimable  if  and  only  if 
(5-7)  L’CX'xrX'X  =  L» 

for  anjr  (X‘X)  “  -  g^X'X)  . 

Proof.  We  show  that  (5-2)  and  (5-7)  are  equivalent.  Clearly  (5*7) 
implies  (5*2);  conversely 

(5.8)  L’(X*X)"X‘X  =  K'X(X'X) "X’X  =  KfX  =  L*  , 

since  X(X'X)"X*X  =  X  [cf.  Pringle  and  Rayner  (1971),  p.  26]. 

Q.E.D. 

We  may  use  (5*7)  to  computationally  verify  estimability  as  follows. 
Substituting  (U.8)  and  (4.9)  into  (5*7)>  with  X*  =  (X'X^X*  gives 
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[l  R-1^ 

(5-9)  L'TT  ~  ~  ~  Tf*  =  L* 

~  -  0  0 

V-  -  J 

Substituting  (5*5)  into  (5*9)  yields 

(5.10)  Lj_R'1S  =  L£ 

To  verify  (5.1C),  therefore,  we  solve  RW  =  S  for  W  ,  say,  which  equals 
R  ,  with  R  upper  triangular.  We  then  examine  L*W  -  and  if  close 
enough  to  0  conclude  L'y  estimable. 

For  the  remainder  of  this  section  we  will  assume  L *7  estimable. 

From  (4.3), 

(5.11)  L'r  =  L'(X'X)~X'y  =  L'X*y  , 

where  X*  =  (X*X)~X*  =  g-,^(X^  ,  cf.  (4.9).  Thus 

B-1  «A 

(5.12)  L'x  =  L'Tt  ~  ~  P*y  =  L’ RXz  , 

~  ~  ~~  0  0/~~  ~  ~ 

V  ~  -7 

A 

using  (4.12)  and  (5.5).  We  compute  L ‘y  ,  therefore,  by  solving  Rw  = 
for  w  ,  say,  which  equals  r’^z^  ,  with  R  upper  triangular.  We  then 
premultiply  by  L£  which  contains  the  r  columns  of  L*  corresponding 
to  the  r  linearly  independert  columns  of  X  which  yielded  R  .  We  note 
that  L*7  is  uniquely  determined  by  (5*1-1)  for  any  (X'X)  =  g^(XfX)  . 

To  see  this,  set  L»  =  K'X  from  (5*2),  so  that  L'(X'X)"X»  =  K*X(X»X)-X'  = 

w  m  M  M  M  M  M  M  M  w 

K,X(X'X)+X*  =  L»(X'X)+X'  ,  since  X(X*X)'X»  is  unique  [cf.  Pringle  and 
Rayner  (1971),  p.  25]. 
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We  define  the  general  linear  hypothesis 


(5.1 l')  L  '7  =  0 

as  testable  whenever  L*7  is  estimable.  The  numerator  of  the  usual  F-test 
for  testable  (5*15)  is  then,  cf.  (2.2), 

(5  *1*0  7,L[L*(X,X)_LrL»7  =  Sh  . 

To  see  that  (5*14)  is  invariant  over  choices  of  (X’X)  ,  notice  that 

L'(X*X)'L  =  K'X(X,X)"X*K  =  K*X(X'X)+X,K  =  L’(X*X)+L  from  (5*2).  Moreover, 
(5.14)  is  also  invariant  over  choices  of  [L'(X'X)  L]  ;  writing 
X*  =  (X’X^X1  we  find  that  (5*1*0  may  be  written 

(5-15)  y'(X*),L[L'X*(X*),L]_L*X*y  =  Sh  , 

using  (5-7)  and  (5. 11).  Sh  is  uniquely  defined  since  for  any  A  , 

A(A’A)  A'  is  unique  [cf.  Pringle  and  Rayner  (1971),  p.  25]. 

To  compute  we  see  from  (4.9)  and  (5*11)  that  (5*15)  may  be  written 

(5.16)  sh  =  *i(H‘1),i1[i[R-V1)*yLiB-135l  • 

We  obtain  an  orthogonal  triangular  decomposition  of 
-1  C\ 

(5*17)  G  =  (R  X),L1  =  Q  ~  ~  ITT*  , 

say,  where  B  is  upper  triangular  txt  ,  with  t  =  r(L)  =  r(L1)  by  (5. 10). 
The  orthogonal  matrix  Q,  is  the  product  of  t  Householder  transformations, 
while  the  permutation  matrix  TTg  rearranges  the  columns  of  ,  r  x  s  , 
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1 


to  make  the  first  t  linearly  independent.  Substituting  (5-17)  into 
(5 .16)  yields 


(5.18) 

* 

S,  =  z|  GG  z, 

h  _  J _ 

y 

where 

2  =  g123^  is 

given  by 

f  -1 

\ 

(5.19) 

*  /B 

°\ 

“  Q'  • 

~  ~  l  2 

Zj 

We  partition  Q  =  ,  where  is  rxt  and  r  x  (r-t)  . 

(if  t  =  r  ,  Q.  -  Q  . ]  Then  (5*18)  reduces  to 

(5-20)  Sh  =  , 

as  at  (2.5).  We  compute  (5*20)  by  applying  the  t  Householder  transfor¬ 
mations  of  Q  in  (5-17)  to  simultaneously  with  G  and  then  summing 

the  squares  of  the  first  t  components  of  the  transformed  z^  . 

If  we  test  the  hypothesis 

(5.21)  L »y  =  m 

and  L'  is  s  x  q  with  row  rank  t  <  s  then  m  must  satisfy  the  same 
s-t  restrictions  that  apply  to  the  rows  of  L'  ,  i.e.,  (5*21)  must  be 
consistent .  Then  the  numerator  cum  of  squares  is  uniquely  given  by 

(5.22)  (y’L-mOll-UX’Xj'Lj'lLv -m)  =  Sh  ; 

following  (5*15)  and  (5*16)  we  see  that 
(5.25)  L’(X»X)"L  =  L|R_1  (R”1)*!^  =  G'G 
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for  which  we  want  a  g^-inverse.  We  use 
LEMMA  $.1.  If  A*  =  g^CA)  ,  then 

(5.24)  A*(A*) «  =  6i2(A»A)  . 

Proof.  From  (4.4),  (4.5)  and  (4.6)  we  have 

(5.25)  A  A*  A  =  A  ,  A*  AA*  =  A*  ,  AA*  =  (A*)  'A'  . 

Hence  A* (A*)  *A'A  =  A*  A  A*A  =  A*A  .  Thus  A'A[A*(A*) 'A^]  =  A' AA*A  =  A'A 
and  [A  (A  )  'A’A]A  (A  ) '  =  A  A  A  (A  )'  =  A  (A  )  »  . 

Q.E.D. 

From  Lemma  5*-  we  obtain 

(5.26)  G*(G*)’  =  [L'(X»X)“L]“ 

from  (5.19),  where  we  partition  tt2  =  (tT2^,TT22)  >  with  E21  *  S  x  t  ’ 
identifying  t  linearly  independent  columns  of  L^  ,  r  x  s  .  Hence 

(5.27)  sh =  ^’5;'™,)I!2ij?"1(!"1),E2i(i,!'!!)  * 

First  L'7  -m  is  computed  and  rearranged  to  form  TT2^(L*7  -m)  =  h  ,  say. 
Then  h  =  B'k  is  solved  for  ,  where  B'  is  lower  triangular.  Finally 
is  found  as  the  sura  of  squares  of  the  components  in 
k  =  (B_1)*h  =  (B-1)  »lg1(L,7  -  m)  • 

The  relevant  F-test  for  the  hypotheses  (5*13)  or  (5*21)  is  then 


computed  as 


(5-28) 


cf.  (2.15),  with  the  critical  region  formed  by  values  of  (5*28)  exceeding 
the  corresponding  tabulated  value  of  F  with  t  and  n-r  degrees  of 
freedom. 

The  above  procedures  simplify  slightly  when  the  contrast  matrix  L’  , 
sxq  ,  has  full  rank  s  <  r  =  r(X)  .  In  that  case  (5.23)  becomes  non¬ 
singular  and  the  results  of  Lemma  5-1  arc  not  needed.  We  use 

LEMMA  5.2-  When  L*7  is  estimable, 

(5.29)  r[L‘(X*X)"L]  =  r(L)  , 


where  r(-)  denotes  rank. 


Proof.  Using  (5-7),  r(L)  =  r[L* (X»X) _X’X]  <  r[L»(X*X)'X» ]  = 
r[L,(X,X)"X,X{(X*X)~},L]  =  r[L*(X'X)"L]  <  r(L)  . 

Q.E.D. 

When  L*  ,  sxq  ,  has  full  row  rank  s  <  r  the  decomposition  (5*17) 
becomes 


(5.30) 


say,  where  tT21  is  now  s  x  s  and  may  equal  Ig  (no  column  interchanges) . 
Formula  (5.27)  applies  with  essentially  no  change. 

We  defer  discussion  of  updating  techniques  for  the  less  than  full  rank 
case  and  extensions  to  multivariate  models  to  a  further  paper.  A  computer 
program  in  Fortran  IV  for  the  IBM  360  is  being  developed  for  the  procedures 
discussed  in  this  paper. 
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